#include <boost/program_options.hpp>

#include <fstream> 
//#include <life/lifepoly/jacobi.hpp>
//#include <life/lifepoly/legendre.hpp>

#include <life/lifepoly/im.hpp>
#include <life/lifepoly/polynomialset.hpp>

#include <math.h>
#include <sstream>
using namespace Life;



template<uint N>
void run() {

//   const uint N = 3;//(const uint) it;
    //const uint N = (const uint)NN;// vm["N"].as<uint>();
 
    detail::OrthonormalPolynomialSet<1,N,Scalar, double,Simplex> d;

    IM<1,20> im;
    ublas::matrix<double> res( d.evaluate( im.points() ) );
    //std::cout<<im.points();
    ublas::vector<double> Pk(N+1);

    for (uint k=0;k<=N;++k) {
      Pk(k)=0;
      for( uint q = 0; q < im.points().size2(); ++ q) {
        Pk(k) += im.weight(q)*res(k,q)*exp(-im.points()(0,q));
      }
    }

    //---------------------------------------------------------------------//

    ublas::vector<double> evalPquad(im.points().size2());
    for (uint q=0;q<im.points().size2();++q) {
      evalPquad(q)=0;
      for (uint k=0;k<=N;++k) {
        evalPquad(q) += Pk(k)*res(k,q);
      }
    }

    double Erreur=0;
    for( uint q = 0; q < im.points().size2(); ++ q) {
        Erreur += im.weight(q)*( exp(-im.points()(0,q)) -evalPquad(q))*( exp(-im.points()(0,q)) -evalPquad(q));
      }
    Erreur=math::sqrt(Erreur);

    std::ofstream fichier1;
  
    fichier1.open("Erreur.txt", std::ios::app); //on ouvrre le fichier en ecriture

    fichier1 << N <<" "<<Erreur<<"\n"; 
    

    fichier1.close();


    //---------------------------------------------------------------------//

    uint nbPts=10;
    ublas::matrix<double> discretisation(1,nbPts+1);
    for( uint q=0;q<=nbPts;++q) { discretisation(0,q)=-1+q*2./nbPts; }

    ublas::matrix<double> evalLk( d.evaluate( discretisation ) );
    //std::cout<<evalLk;

    ublas::vector<double> evalPq(nbPts+1);
    for (uint q=0;q<=nbPts;++q) {
      evalPq(q)=0;
      for (uint k=0;k<=N;++k) {
        evalPq(q) += Pk(k)*evalLk(k,q);
      }
    }

    std::ostringstream Nstring;Nstring << N;
    std::string name("part2_N"+Nstring.str() );
    
    std::ofstream fichier;
  
    fichier.open(name.c_str(), std::ios::out); //on ouvrre le fichier en ecriture

    for (uint q=0;q<=nbPts;++q) {
      fichier << discretisation(0,q) <<" "<<exp(-discretisation(0,q))<<" "<<evalPq(q)<<"\n"; 
    }

    fichier.close();

  }


int main(int argc, char** argv )
{

run<1>();
run<2>();
run<3>();
run<4>();
run<5>();
run<6>();
run<7>();
run<8>();
run<9>();

return 0;

}

